home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpklib.zip / SPOCO.FOR < prev    next >
Text File  |  1984-01-06  |  6KB  |  195 lines
  1.       SUBROUTINE SPOCO(A,LDA,N,RCOND,Z,INFO)
  2.       INTEGER LDA,N,INFO
  3.       REAL A(LDA,1),Z(1)
  4.       REAL RCOND
  5. C
  6. C     SPOCO FACTORS A REAL SYMMETRIC POSITIVE DEFINITE
  7. C     MATRIX AND ESTIMATES THE CONDITION OF THE MATRIX.
  8. C
  9. C     IF  RCOND  IS NOT NEEDED, SPOFA IS SLIGHTLY FASTER.
  10. C     TO SOLVE  A*X = B , FOLLOW SPOCO BY SPOSL.
  11. C     TO COMPUTE  INVERSE(A)*C , FOLLOW SPOCO BY SPOSL.
  12. C     TO COMPUTE  DETERMINANT(A) , FOLLOW SPOCO BY SPODI.
  13. C     TO COMPUTE  INVERSE(A) , FOLLOW SPOCO BY SPODI.
  14. C
  15. C     ON ENTRY
  16. C
  17. C        A       REAL(LDA, N)
  18. C                THE SYMMETRIC MATRIX TO BE FACTORED.  ONLY THE
  19. C                DIAGONAL AND UPPER TRIANGLE ARE USED.
  20. C
  21. C        LDA     INTEGER
  22. C                THE LEADING DIMENSION OF THE ARRAY  A .
  23. C
  24. C        N       INTEGER
  25. C                THE ORDER OF THE MATRIX  A .
  26. C
  27. C     ON RETURN
  28. C
  29. C        A       AN UPPER TRIANGULAR MATRIX  R  SO THAT  A = TRANS(R)*R
  30. C                WHERE  TRANS(R)  IS THE TRANSPOSE.
  31. C                THE STRICT LOWER TRIANGLE IS UNALTERED.
  32. C                IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
  33. C
  34. C        RCOND   REAL
  35. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .
  36. C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS
  37. C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE
  38. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
  39. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
  40. C                           1.0 + RCOND .EQ. 1.0
  41. C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING
  42. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
  43. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
  44. C                UNDERFLOWS.  IF INFO .NE. 0 , RCOND IS UNCHANGED.
  45. C
  46. C        Z       REAL(N)
  47. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
  48. C                IF  A  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS
  49. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
  50. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
  51. C                IF  INFO .NE. 0 , Z  IS UNCHANGED.
  52. C
  53. C        INFO    INTEGER
  54. C                = 0  FOR NORMAL RETURN.
  55. C                = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
  56. C                     OF ORDER  K  IS NOT POSITIVE DEFINITE.
  57. C
  58. C     LINPACK.  THIS VERSION DATED 08/14/78 .
  59. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  60. C
  61. C     SUBROUTINES AND FUNCTIONS
  62. C
  63. C     LINPACK SPOFA
  64. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  65. C     FORTRAN ABS,AMAX1,REAL,SIGN
  66. C
  67. C     INTERNAL VARIABLES
  68. C
  69.       REAL SDOT,EK,T,WK,WKM
  70.       REAL ANORM,S,SASUM,SM,YNORM
  71.       INTEGER I,J,JM1,K,KB,KP1
  72. C
  73. C
  74. C     FIND NORM OF A USING ONLY UPPER HALF
  75. C
  76.       DO 30 J = 1, N
  77.          Z(J) = SASUM(J,A(1,J),1)
  78.          JM1 = J - 1
  79.          IF (JM1 .LT. 1) GO TO 20
  80.          DO 10 I = 1, JM1
  81.             Z(I) = Z(I) + ABS(A(I,J))
  82.    10    CONTINUE
  83.    20    CONTINUE
  84.    30 CONTINUE
  85.       ANORM = 0.0E0
  86.       DO 40 J = 1, N
  87.          ANORM = AMAX1(ANORM,Z(J))
  88.    40 CONTINUE
  89. C
  90. C     FACTOR
  91. C
  92.       CALL SPOFA(A,LDA,N,INFO)
  93.       IF (INFO .NE. 0) GO TO 180
  94. C
  95. C        RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
  96. C        ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
  97. C        THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
  98. C        GROWTH IN THE ELEMENTS OF W  WHERE  TRANS(R)*W = E .
  99. C        THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
  100. C
  101. C        SOLVE TRANS(R)*W = E
  102. C
  103.          EK = 1.0E0
  104.          DO 50 J = 1, N
  105.             Z(J) = 0.0E0
  106.    50    CONTINUE
  107.          DO 110 K = 1, N
  108.             IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K))
  109.             IF (ABS(EK-Z(K)) .LE. A(K,K)) GO TO 60
  110.                S = A(K,K)/ABS(EK-Z(K))
  111.                CALL SSCAL(N,S,Z,1)
  112.                EK = S*EK
  113.    60       CONTINUE
  114.             WK = EK - Z(K)
  115.             WKM = -EK - Z(K)
  116.             S = ABS(WK)
  117.             SM = ABS(WKM)
  118.             WK = WK/A(K,K)
  119.             WKM = WKM/A(K,K)
  120.             KP1 = K + 1
  121.             IF (KP1 .GT. N) GO TO 100
  122.                DO 70 J = KP1, N
  123.                   SM = SM + ABS(Z(J)+WKM*A(K,J))
  124.                   Z(J) = Z(J) + WK*A(K,J)
  125.                   S = S + ABS(Z(J))
  126.    70          CONTINUE
  127.                IF (S .GE. SM) GO TO 90
  128.                   T = WKM - WK
  129.                   WK = WKM
  130.                   DO 80 J = KP1, N
  131.                      Z(J) = Z(J) + T*A(K,J)
  132.    80             CONTINUE
  133.    90          CONTINUE
  134.   100       CONTINUE
  135.             Z(K) = WK
  136.   110    CONTINUE
  137.          S = 1.0E0/SASUM(N,Z,1)
  138.          CALL SSCAL(N,S,Z,1)
  139. C
  140. C        SOLVE R*Y = W
  141. C
  142.          DO 130 KB = 1, N
  143.             K = N + 1 - KB
  144.             IF (ABS(Z(K)) .LE. A(K,K)) GO TO 120
  145.                S = A(K,K)/ABS(Z(K))
  146.                CALL SSCAL(N,S,Z,1)
  147.   120       CONTINUE
  148.             Z(K) = Z(K)/A(K,K)
  149.             T = -Z(K)
  150.             CALL SAXPY(K-1,T,A(1,K),1,Z(1),1)
  151.   130    CONTINUE
  152.          S = 1.0E0/SASUM(N,Z,1)
  153.          CALL SSCAL(N,S,Z,1)
  154. C
  155.          YNORM = 1.0E0
  156. C
  157. C        SOLVE TRANS(R)*V = Y
  158. C
  159.          DO 150 K = 1, N
  160.             Z(K) = Z(K) - SDOT(K-1,A(1,K),1,Z(1),1)
  161.             IF (ABS(Z(K)) .LE. A(K,K)) GO TO 140
  162.                S = A(K,K)/ABS(Z(K))
  163.                CALL SSCAL(N,S,Z,1)
  164.                YNORM = S*YNORM
  165.   140       CONTINUE
  166.             Z(K) = Z(K)/A(K,K)
  167.   150    CONTINUE
  168.          S = 1.0E0/SASUM(N,Z,1)
  169.          CALL SSCAL(N,S,Z,1)
  170.          YNORM = S*YNORM
  171. C
  172. C        SOLVE R*Z = V
  173. C
  174.          DO 170 KB = 1, N
  175.             K = N + 1 - KB
  176.             IF (ABS(Z(K)) .LE. A(K,K)) GO TO 160
  177.                S = A(K,K)/ABS(Z(K))
  178.                CALL SSCAL(N,S,Z,1)
  179.                YNORM = S*YNORM
  180.   160       CONTINUE
  181.             Z(K) = Z(K)/A(K,K)
  182.             T = -Z(K)
  183.             CALL SAXPY(K-1,T,A(1,K),1,Z(1),1)
  184.   170    CONTINUE
  185. C        MAKE ZNORM = 1.0
  186.          S = 1.0E0/SASUM(N,Z,1)
  187.          CALL SSCAL(N,S,Z,1)
  188.          YNORM = S*YNORM
  189. C
  190.          IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
  191.          IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
  192.   180 CONTINUE
  193.       RETURN
  194.       END
  195.